%--------------------------------------------------------------------------
% Computes the basic state solution of the problem
%
%             v_t = v_{xx} 
% subject to:
%
% - no flux bc at bottom
% - evaporation at top
% - v = 0 at t = 0
%
% note that here, v = c - beta, where c is the concentration and beta is
% the initial condition
%
%--------------------------------------------------------------------------

function soln = base(p)

if nargin == 0
    p = params;
end

eta = linspace(0, 1, p.N);


X = [zeros(p.N, 1); 1];

% use the steady state as the initial condition
% v = -p.delta * p.beta * (1 - p.beta) * (eta.^2 / 2 - 1/6)';
% X = [v; 1];

opts = odeset('AbsTol',1e-8, 'RelTol',1e-8);
soln = ode15s(@(t, X) system(t, X, p), [0, p.Tmax], X, opts);


function F = system(t, X, p)

v = X(1:p.N);
h = X(end);

eta = linspace(0, 1, p.N)';

dvdn = comp_base_deriv(v, h, p);

d2vdn = [
    (2 * v(2) - 2 * v(1)) / p.dz^2;
    (v(1:end-2) - 2 * v(2:end-1) + v(3:end)) / p.dz^2;
    (2 * v(p.N-1) - 2 * v(p.N) + 2 * p.dz * dvdn(p.N)) / p.dz^2
    ];

dh = -p.delta * (p.beta + v(p.N));

% dvdn(1) = 0;
% d2vdn(1) = 0;

F = [eta * dh .* dvdn / h + d2vdn / h^2; dh];